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ABSTRACT 

In this paper, we present a model for the long-term evolution of the merger of two unequal mass C/O 
white dwarfs (WDs). After the dynamical phase of the merger, magnetic stresses rapidly redistribute 
angular momentum, leading to nearly solid-body rotation on a viscous timescale of 10~ 4 — 1 yr, long 
before significant cooling can occur. Due to heating during the dynamical and viscous phases, the less 
massive WD is transformed into a hot, slowly rotating, and radially extended envelope supported by 
thermal pressure. 

Following the viscous phase of evolution, the maximum temperature near the envelope base may 
already be high enough to begin off-center convective carbon-burning. If not, Kelvin-Helmholtz con- 
traction of the inner region of the envelope on a thermal timescale of 10 3 — 10 4 yr compresses the 
base of the envelope, again yielding off-center burning. As a result, the long-term evolution of the 
merger remnant is similar to that seen in previous calculations: the burning shell diffuses inwards 
over ~ 10 4 yr, eventually yielding a high- mass O/Ne WD or a collapse to a neutron star, rather than 
a Type la supernova. During the cooling and shell-burning phases, the merger remnant radiates near 
the Eddington limit. Given the double WD merger rate of a few per 1000 yr, a few dozen of these 
~ 10 38 erg s" 1 sources should exist in a Milky Way- type galaxy. 

While the end result is similar to that of previous studies, the physical picture and the dynamical 
state of the matter in our model differ from previous work. Furthermore, substantial remaining 
uncertainties related to the convective structure near the photosphere and mass loss during the thermal 
evolution may significantly affect our conclusions. Thus, future work within the context of the physical 
model presented here is required to better address the eventual fate of double WD mergers, including 
those for which one or both of the components is a He WD. 

Subject headings: binaries: close — nuclear reactions, nucleosynthesis, abundances — supernovae: 
general — white dwarfs 



1. INTRODUCTION 



Ty pe la supernovae (SNe la; iHillebrandt fe Niemeverl 
2000) are well known for their role in determining 
the accelerating expansion o f the Universe (jRiess et al.l 
119981: iPerlmutter et al.l I1999T ). While the generally ac- 
cepted model of a thermonuclear runaway in a C/O 
white dwarf (WD) explains many of the defining char- 
acteristics of SNe la, substantial uncertainty remains 
as to the progenitor system. Three progenitor scenar- 
ios are current ly considered: (1) the single degene rate 
scenario (SDS; IWhelan fc Ibenlll973l : INomotol fl98l . in 
which a WD accretes and stably burns H-rich mate- 
rial from a non-degenerate main sequence, sub- giant, 
or red giant donor; (2) the double degenerate scenario, 
in which two C/O WDs undergo a dynamical merger 
(llben fe Tut ukovlll984l: IWebbinkl| 1984 [) : and (3) the dou- 
ble detonat i on scenario (ILivnel I199C : ILivne fe Glasner 
19901 [19911 IWooslev fe Weaverl [1991 ILivne fe Arnett 



1995), in which a detonation in an overlying He-burning 



shell sends a converging shock wave into the C/O core 
that initiates a subsequent detonation near the center. 
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In the SDS and double degenerate scenario, convec- 
tive C-burning is ignited by strongly-screened or py- 
cnonuclear reactions induced by very high densities 
~ 10 10 g cm -3 in the center of a near-Chandra s ekhar 
mass WD. After 3 00 - 1000 yr (|Wooslev et al.l [2001 
iPiro fe Ch ang 2008j), the burning in this convective re- 
gion becomes violent enough to yield a deflagration 
and/or detonation, incinerating the WD. In the dou- 
ble detonation scenario, no core convective burning takes 
place: the shock from a He shell detonation initiates a 
secondary detonation in the C/O core that disrupts the 
WD. While these scenarios cover a wide range of possi- 
bilities, none can adequately explain specific features of 
SNe la and their observed rate. 

The SDS relies on the stable burning of accreted hy- 
drogen in order to g row the WD . How e ver, it ha s been 
known for deca des ([Sienkiewiczl 119801 : iFuiimotol 119821: 
Paczvri ski 198 3|) and recen t ly confirmed ( YoonetaLl 
20041 IShen fc Bildstenl [2007t INomoto et all 120071 ) that 



the range of accretion rates that allow for stable H- 
burning is very narrow. Furthermore, the accumu- 
lated He layer will undergo nova-like flashes that likely 
eject some of the accreted mate rial. As a result , 
population synthesis ca l culations ([Branch et al.l Il995l ; 
lYungelson fc Livid [l99l [20001 IRuiter et al.l I2009D pre- 
dict SN la rates from the SDS channe l that are ~ 10 
times lower than the obseryed rate (IMannucci et al.l 
[2001 [Scannapieco fc Bildstenl [20051 iDildav et al.l 12008ft . 
Additionally, observational limits on the progeni- 
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tor systems (IDi St efano et al. 200 61 IDi Stefand 120101: 

IT 



IGilfanov fc Bogdanl 120101: ILi et al.l [201 lft and on the 
SN la e.jecta's interaction with a donor ( Mariett a et all 
2000HKasenl2010HHavden et al.lfe OlO: Bia nco et al.l2011b 
Nugent et al.H2011t IBloom et al. 120121 ) severely constrain 



the red giant donor channel's contribution to the SN la 
population. 

The double detonation scenario has received more 
attention recently due to the pr ediction of He deto- 
nations in pre- AM CVn systems (IBildsten et al.l 120071: 
Fink et al.ll2007L 1201(1 IShen fc Bildstenl|2009t , IShen et al l 
20101: ISim et allboiOHWooslev fc Kasenl2011l) . However, 
while many features of this scenario are attractive, it is 
not clear how the He-detonation products on the out- 
side of the core will skew observations away from SN 
la-like events. Even a small amount of surface iron- 
group elements created during He-burning, especially Ti 
and Cr, can significantly impact the colors of the SN 
at peak light. The UV line blanketing due to these 
elements causes the predicted spe ctra near peak light 
to be much red der than observed (|Kromer et al.1 1201(1 
ISim et a"ill201lD FI Furthermore, the He detonation will 
eject radioactive isotopes at high velocities, a feature that 
has never been observed in SNe la. 

1.1. An updated model for double white dwarf mergers 

The double degenerate scenario involves the dynam- 
ical coalescence of two C/O WDs. The predicted 
merger rate and delay tim e distribution (IBranch et al.1 
1199a INelemans et alj|200U IRuiter et all 120091) and~oi> 
servations of progenitor sy stems (albeit with large e r- 
ror bars: INelemans etaLl 120051: iGeier et all l2010allbl ; 



iTovmassian et al.l 120101: iRodriguez-Gil et al.l (20101) are 
promisingly consistent with t he observed SN la rate 
and delay time distribution (jMaoz et al.l 120101 120111 : 
IGraur et al.ll201lh . However, previous theoretical stud- 
ies have found that the merger of two WDs in- 
stead triggers an accretion induced collapse, yielding 
a faint explosi on and a neutron star remnant instead 
of a SN la (INomoto fc Ibenl [19851 IKawai et al.1 H987I: 
ISaio fc Nomotdl2004D . 

These studies assumed that the less massive WD 
is disrupted into a Keplerian disk from which the 
more massive WD accretes at a constant rate near 
the gravitational Eddington limit of 10~ 5 M© yr -1 
for ~ 10 5 yr. In contrast, smoothed particle hy- 
drodynamical (SPH) studies of the merging process 
(iBenz et al.lll990l:lRasio fc Shapirdll995HSegretain et all 
1997 1 : iGuerrero et all 120041: lLoren-Aguilar et all 120051 
2009[ lYoon et all 120071: iDan et all 120111: iRaskin et ail 
20111) find that the tidally disrupted WD is instead con- 



verted into material that is supported against gravity by 
a combination of thermal pressure from shock heating 
in regions close to the more massive WD and centrifu- 
gal support in the optically th ick o uter regions. Until 
the w ork of lYoon et all (|2007l ) and Ivan K crkwii k et all 
(2010), these important features of the merged configu- 
ration were largely neglected in studies of the long-term 
evolution of merger remnants. 

One conclusion reached from consideration of this more 

6 Conversely, previous studies found t he colors to be too blue due 
to 56 Ni decays ne ar the photosphere (Hoflich & Khokhlovl 119961 : 
INugent et alt[T997l ). 



detailed merger end-state is that the viscous evolution of 
the differentially rotating layers must be taken into ac- 
count before the secular thermal evolution of the merger 
remnant can be calculated. In particular, if magnetically- 
govern ed processes such as the magnetorotational insta- 
bility dBalbus fc HawlevI 119911: lHawlev fc Balbus) 119911: 



i BalbusI 120031) and the Tayler-Spruit dynamo ( Tavlerl 
ll973HSpruitHl999l 12002ft operate efficiently, angular mo- 
mentum will be redistributed throughout the remnant 
much more rapidly than heat can be radiated away. 
Thus, the envelope consisting of the tidally disrupted 
WD will transfer most of its angular momentum to a 
small amount of mass in the outer layers and approach a 
viscously-heated shear-free configuration long before cool- 
ing significantly. 

In this paper, we present an updated physical model 
for the evolution of unequal mass double WD mergers. 
Dynamical and viscous processes, acting on timescales 
much shorter than the thermal timescale, convert the 
majority of the less massive WD into a hot radially ex- 
tended envelope surrounding the more massive WD. Lit- 
tle mass remains centrifugally-supported following the 
viscous evolution. While the entire remnant has a mass 
essentially equal to the total mass of the initial binary 
and may be larger than the Chandrasekhar mass, the 
hot envelope is extended and does not exert considerable 
pressure on the degenerate core. 

Substantial heat is generated during the dynamical and 
viscous phases of evolution and may lead to convective 
C-burning at the interface between the hot envelope and 
degenerate core during or directly after the viscous evo- 
lution. If instead the maximum temperature is too low 
for immediate convective C-burning prior to significant 
thermal transport, the inner part of the envelope will 
undergo a Kelvin-Helmholtz contraction on a thermal 
timescale of 10 3 — 10 4 yr, which is set by the Eddington 
luminosity of the giant envelope. The hot non-rotating 
envelope radiates near the Eddington limit, cooling and 
compressing the base of the envelope relatively rapidly, 
just as the base of the accreted envelope is compressed in 
studie s assuming disk accretion such as INomoto fc Ibenl 
(1985). This contraction proceeds on a shorter timescale 
than the inward thermal diffusion timescale into the de- 
generate core of ~ 10 6 yr, so that a temperature inversion 
is formed and a convective C-burning shell is born well 
away from the center. 

Thus, it appears unavoidable that double C/O WD 
mergers ignite C-burning well away from the merger rem- 
nant's center. This off-center ignition allows the convec- 
tive layer to expand, and the necessary conditions for a 
violent deflagration and/or detonation are not reached. 
Under the assumption of heat transport by thermal con- 
duction, the burning shell moves slowly inwards at a 
speed of ~ 0.01 cm s _1 , reaches the center in ~ 10 4 yr, 
and eventua lly converts the prev i ously C/O WD into a n 
O/Ne WD (|Timmes et allll99l ISaio fc Nomotol [1991 . 
If the remnant mass is larger than the Chandrasekhar 
mass, envelope material can continue to contract and 
compress the O/Ne WD until it collapses into a neu- 
tron star, accompanied by a weak explosion. While 
such events will not yield SNe la, they may give rise to 
fast transients and contribute to r-proce ss enrichment, 
and are certainly of astroph y sical interest (iDessart et all 
I2006t iMetzger et all l2009a.b; D arbha et al.ll2010ft . 
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Although this chain of events represents an updated 
description of double WD merger evolution, the final out- 
come is similar to that found by previous studies assum- 
ing disk accretion. The onset of convective C-burning 
in a shell following the viscous phase of evolution, cither 
prior to or after a period of thermal redistribution, seems 
inevitable. However, the physical picture leading to this 
evolution is markedly different from previous work. In 
our model, there is no phase of accretion from a Keple- 
rian disk onto a degenerate object; instead, the long-term 
thermal evolution of the merger remnant is determined by 
the conditions at the base of the hot non-rotating envelope 
and the rate at which it redistributes heat and cools. 

In this paper, we describe our updated model for the 
evolution of a double C/O WD merger, focusing sepa- 
rately on the dynamical (jJ5J, viscous (fj3|), and thermal 
O phases. In fj5j we compare to previous work. We 
conclude in fj6] and speculate on future work and poten- 
tial observables. 

2. DYNAMICAL EVOLUTION: 10 2 - 10 3 s 

Close binaries consisting of tw o C/O WDs are forme d 
via two common envelope phases (jlben fe T utukov 1984): 
one while the initially more massive star is on its asymp- 
totic giant branch (AGB), and another during the second 
star's AGB or at the tip of its red giant branch. In the 
latter case, the second star lives for ~ 100 Myr as a core 
He-burning s dO or sdB star before becoming a C/O WD 
(|Heberll2009j ). The orbital separations of double WD bi- 
naries shrink via gravitational wave radiation, bringing 
WD binaries with initial orbital periods < 10 hr into 
mass transfer contact within a Hubble time. These bina- 
ries are being discovered in growing numbers by studies 
such as the SPY and the SPSS-based SWARMS surveys 
(I York et al.M2000t INelemans et all 120051: iBadenes et al.l 
l2009HMullallv et al.H2009HGeier et al.ll2010afl . 

Upon contact, the less massive WD overflows its Roche 
lobe and begins transferring matter onto the more mas- 
sive WD. The dynamical stability of this mass transfer 
depends on the mass ratio and on how efficiently angular 
mo mentum is returned to the orbit during mass trans- 
fer (jMarsh et al.H200"4T ). For typical assumptions, mass 
transfer between two C/O WDs is dynamically unstable, 
even for the lowest expected mass ratio of 0.5/1.2. 

The unstable mass transfer is driven by the less mas- 
sive WD's hydrodynamical response to mass loss, and 
thus the dynamical evolution is characterized by that 
WD's sound crossing time or, nearly equivalently, the 
binary's orbital period: 30 — 300 s. The subsequent evo- 
lution of the resulting merger remnant is dominated by 
the transport of angular momentum and entropy. Sim- 
ple estimates of the importance of magnetic stresses lead 
to viscous evolution within 10 4 — 10 8 s. This timescale 
is much longer than the dynamical timescale, but much 
shorter than the remnant's thermal timescale, 10 3 — 10 4 
yr. Thus, we separate discussion of the merger and the 
subsequent evolution of the remnant into three distinct 
phases based on this hierarchy of timescales. 

The dynamical merging p rocess has been stud- 
ied i n detail (IBenz et al.l Il990t IRasio fc Shapiro] 
1995t iGuerrero et al.l 12004 Loren-Aguilar et al.l 120091: 
Dan et al. 120111 : Ras kin et al.l l 20111) . so we only empha- 
size a few results. Once mass transfer begins, material 
leaves the less massive WD's Roche nozzle and directly 




M enc (r cyl )/M e from center of 0.9 M WD 



Figure 1. Angular velocity for material within 10 s cm of the equa- 
torial midplane (solid line), Keplerian angular velocity (dashed 
line), and their ratio (dotted line) vs. cylindrically enclosed mass 
from the center of the 0.9 Mq WD. Data are taken at the end of 
an SPH simulation of a merging 0.6 + 0.9 Mq WD binary, 360 s 
6 initial or bital periods) after the 0.6 Mq WD has been disrupted 
Dan et al.ll20Tll) . 

impacts the surface of t he more massive WD (IWebbinki 
ll984HMarsh et al.ll2004T) . The resulting shock efficiently 
converts the gravitational and kinetic energy into ther- 
mal energy, so that this material is nearly virializcd. As 
tidal disruption progresses in earnest, the less massive 
WD becomes unbound from itself. Because much of 
this material has not been shocked, its entropy stays 
relatively low, and it remains supported against gravity 
by centrifugal forces as it orbits the more massive WD. 

As a fiducial example, we show results from an SPH 
simulation of a me rging 0.6 + 0.9 M P) double WD bi- 
nary performed bv iDan et alJ (|2011h . which is qualita- 
tively representative of other SPH calculations of unequal 
mass merge rs. The WDs are as sumed to be rotating syn- 
chronously (|Fuller fc Laill2011l) . The solid line in Figure 
Q] shows the angular velocity, f2, of material within 10 8 
cm of the equatorial midplane, taken with respect to the 
center of the 0.9 Mq WD, vs. total cylindrically en- 
closed mass 360 s after the 0.6 Mq WD has been dis- 
rupted. Also shown is the Keplerian angular velocity, 
= (GMenc/r^yj) 1 / 2 (dashed line), where M enc is the 
mass enclosed within a sphere of radius r cy \. The dotted 
line shows their ratio. Following the merger, centrifugal 
forces are the dominant support against gravity for the 
outer 0.4 Mq of material. Note that the outer 0.1 Mq 
exists as a spiral tidal tail and thus does not appear fully 
Keplerian. 

Figure [2] shows the enclosed mass vs. spherical radius 
from the center of the 0.9 Mq WD. The solid line shows 
the results from the beginning of the simulation, just 
before mass transfer begins. Mass transfer continues at 
increasing levels for 1840 s (30 initial orbital periods) 
until the less massive WD is disrupted. The dashed line 
shows the enclosed mass vs. radius an additional 360 s 
after the 0.6 Mq WD has become unbound. 

The enclosed mass vs. radius for most of the inner 
0.9 Mq of the remnant remains largely undisturbed be- 
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?"sph [cm] from center of 0.9 M Q WD 

Figure 2. Spherically enclos ed mass vs. spher ical radius from the 
center of the 0.9 Mq WD in IDan et all puTfl Vs SPH simulation. 
The solid line shows results from the beginning of the simulation, 
just before mass transfer begins. The dashed line shows results 
360 s after the 0.6 Mq WD has been disrupted. Also shown is the 
zero- temperature mass-radius relation for C/O WDs (dotted line). 

cause most of the tidally disrupted WD remains sup- 
ported against gravity by rotation, with a smaller contri- 
bution from thermal pressure. However, a small amount 
of the disrupted WD does exert some pressure on the 
0.9 Mq core. The dotted line shows the zero-temperature 
mass-radius relation for C/O WDs, which intersects the 
"before" (solid) and "after" (dashed) lines at masses that 
differ by 0.05 M Q . 

Some recent studies have focused on the dynamical 
mergers of binaries with equal or nearly equal mass 
components. In contrast to the SPH results we use 
throughout this paper, these studies have found that 
both co mponents, instead o f just one, are tidally dis- 
rupte d (jLoren-Aguilar et"al~l 120091 : Ivan Kerkwijk et aTl 
I2010D or that the resulting merger is more violent and 
reach es higher temperatures th an in the unequal mass 
case (|Pakmor et al.|[2~010L l201l[) . However, it is possible 
that these differences are artifacts of the approximate ini- 
tial conditions these particular equal mass studies use. 
Namely, these simulations begin with the binary com- 
ponents close enough that strong mass transfer immedi- 
ately sets i n onc e the calculation is begun. In contrast, 
IDan e t al. (2011) emphasize the importance of beginning 
such simulations at larger orbital separations and instead 
find tidal disruption at a much larger radius with corre- 
spondingly less violence. Given the current uncertainties, 
we proceed under the assumption that only one WD is 
tidally disrupted during the merger. We further assume 
that the extreme conditions n ecessary for a detonation 
are not met during the merger (jPakmor et al.ll2lMl20l"l"l : 
IGuillochon e till 120101) . 

3. VISCOUS EVOLUTION: 10 4 - 10 8 s 

Once the merger remnant has evolved for several dy- 
namical timescales, it will reach quasi-hydrostatic equi- 
librium. As Figure [1] shows, considerable differential ro- 
tation exists in the disrupted WD's material, more than 
half of which is nearly Keplerian. This shearing in the 



presence of viscous mechanisms naturally leads to out- 
ward angular momentum transport. Quantitative re- 
sults of this viscous evolution require a multi-dimensional 
magnetohydrodynamic calculation; suc h two- and three- 
dime nsional simulations with ZEUS (jStone fe Normanl 
1992) are currently underway. However, some general 
points can be made. 

Application of an a-viscosity prescription 
(Shak ura fc Sunvaevj I1973D in the geometrically 
thick near-Keplerian materia l, e.g. due to th e mag- 
netorotational instabi l ity dBalbus fc Hawleyl 119911 : 
lHawlev fc Balbusl Il99lt iBalbusI I2003D . yields a viscous 
timescale of 




where r cy i is the cylindrical radius, and a typical pres- 
sure scale height in this rotating material is h ~ 0.1r cy i. 
Magnetic stresses may also yield efficient angular mo- 
mentum transport in the non-rotationally-supported ma- 
terial that has been spun up during the merger (e.g., the 
material between 0.9 — 1.0 M Q in Fig. [T]). While the mag- 
netorotational instability cannot operate in this regime 
where dCl/dr > 0, stresses due to the winding up of a 
seed magnetic field by the Tayler-Spruit dynamo may b e 
very effective (|Tavlerft 19731: lSpruitlll999l l2002t lPircl2008f ) . 
Since this material is stably stratified by its thermal pr o- 
file, the appropriate effective diffusivity is lSpruitl (|2002Tj 's 
case 1, 

- 3 x 10 io ( r ^ \ ( 9. y /2 

J g y'V 18 ^ (2) 
V3x 10 5 cm 2 /sy \ N )' K ' 

The square of the Brunt- Vaisala, frequency in the absence 
of compositional gradients is 

where the gravitational acceleration is g, the logarithmic 
derivative of the pressure is Xx = d In P/d In x with the 
other intrinsic variables held constant, and the adiabatic 
temperature gradient is V a d = 91nT/(31nP at constant 
entropy. The thermal diffusivity in equation (j2J is given 

by 

16crT 3 

3np 2 cp 
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0.01 cm 2 /g\ fl.5k/Um p 



cp 



(4) 



where the opacity is K, and the specific heat at constant 
pressure, cp, is due primarily to ions in this mostly de- 
generate material. This yields a viscous timescale of 

,.2 



^visc.TS ' 



^5 x 10 7 s 



r cy i 



n 



3 x 10 s cm/ V0-3S- 1 
3 x 10 10 cm 2 /s\ /0.02 s" 1 



(5) 



where the shear rate in these layers is q = dVl/ d\nr cy \. 

As long as a <C 1, both the a-viscosity and Tayler- 
Spruit dynamo processes yield viscous timescales that 
are 10 4 — 10 8 s, so that our assumption that the merger 
remnant enters into hydrostatic equilibrium prior to its 
viscous evolution is well-founded. Furthermore, as we 
show in the next section, the viscous timescale is orders 
of magnitude shorter than the thermal timescale. Thus, 
the viscous evolution occurs with no significant heat loss. 
Because the material is optically thick and neutrino cool- 
ing is not significant on such short timescales, the heat 
generated dynamically during the merger and through 
viscous shear heating will remain trapped in the enve- 
lope until it can diffuse out on the much longer thermal 
timescale. 

As angular momentum is transported outwards to 
larger radii, the merger remnant evolves towards an equi- 
librium state of shear-free solid-body rotation. During 
the progression towards this equilibrium, the local spe- 
cific heating rate due to shear, e V isc = vq 2 , both in- 
creases the specific internal energy, u, and does expan- 
sion work by spreading material to larger radii. While 
the precise end state of the remnant is a function of 
the viscous mechanism and the multi-dimensional time- 
dependent evolution, local conservation of energy im- 
plies that the originally near-Keplerian material is con- 
verted to nearly virial material with total specific en- 
ergy u — GM cnc /r ~ (r cy if2o) 2 /2, where fin is the ini- 
tial angular velocity after the dynamical phase of evo- 
lution. In effect, during the course of viscous evolution, 
the centrifugally-supported material is converted into a 
hot envelope surrounding the degenerate core. By the 
time significant thermal transport begins, most of the 
less massive WD has been converted into a thermally- 
supported and extended envelope with negligible rota- 
tion. 

As an example, we present a one-dimensional a- 
vi scosity hyd r odyna mic simulation of the outer 0.6 M Q 
of iDan et afl (j2011f l's 0.6 + 0.9 M Q merger. Details of 
this viscous calculation are presented in the Appendix. 
The simulation begins 360 s after the 0.6 M© WD has 
been disrupted. We assume the inner 0.9 Mq remains 
unchanged during the viscous evolution. We evolve the 
axisymmetric height-integrated Navier-Stokes equations 
under the assumption that the only entropy evolution 
term is due to shear heating because of the inefficiency 
of radiative cooling. We follow the cylindrical radius, 
r cy i, radial and angular velocities, v r and f2, height- 
integrated surface density, £ = pdz, and height- 
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Figure 3. Measures of the angular velocity vs. cylindrically en- 
closed mass at different times during the ID viscous simulation. 
Solid lines show the angular velocity, and dotted lines show the 
ratio of at 0, 10 3 , and 10 4 s after the beginning of the sim- 

ulation. 

integrated pressure multiplied by a factor of the adia- 
batic exponent 7, II = 7 J_ Pdz, for 200 equal mass 
Lagrangian cylindrical shells. 

We assume a — 10~ 2 throughout the material, even 
though the physical source of viscosity changes from the 
inner regions to the outer. While this could alter the pre- 
cise values at the end of the simulation, the qualitative 
results will be similar as long as the viscous timescale is 
intermediate to the dynamical and thermal timescales. 
A 7 = 5/3 equation of state is assumed to hold through- 
out the disk, which becomes questionable when radiation 
pressure becomes important at later times. However, 
simulations with smaller values of 7 have qualitatively 
similar outcomes. 

Figure [3] shows the resulting angular velocity vs. the 
cylindrically enclosed mass 0, 10 3 , and 10 4 s after the 
beginning of our simulation (solid lines). As expected 
from equation (fTJ), the rotational profile evolves towards a 
shear-free configuration within a viscous timescale ~ 10 4 
s. Also plotted are the ratios of SI/VLk at these different 
times (dotted lines), showing that the material originally 
from the tidally disrupted WD evolves from being mostly 
rotationally-supported to mostly pressure-supported. 

Figure 3] shows the specific angular momentum pro- 
file vs. cylindrical radius at the same snapshots in time 
as in Figure [3] Markers denote cylindrically enclosed 
masses as labeled. It is clear that as angular momen- 
tum is transferred outwards, the outer regions > 1.2 Mq 
evolve by spreading. Figure [5] shows the ratio of II/E 
vs. cylindrically enclosed mass at the same snapshots in 
time as in Figure [3] This ratio is a proxy for the spe- 
cific internal energy at the equator; calculation of the 
actual equatorial specific energy requires knowledge of 
the remnant's vertical structure, which is not followed 
in the height-integrated approximation. As Figure [5] 
qualitatively demonstrates, the internal energy between 
1.0 — 1.3 Mq increases significantly as rotational energy 
is converted into heat. Because of this thermal pressure 
support, material within 1.3 Mq does not move inwards 
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Figure 4. Specific angular momentum vs. cylindrically enclosed 
radius at the same snapshots in time as in Figure[3] Markers denote 
cylindrically enclosed masses of 1.0 — 1.4 Mq separated by 0.1 M , 
as labeled for the 10 4 s snapshot. 
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Figure 5. A proxy for the equatorial specific internal energy, 
II/S, vs. cylindrically enclosed mass at the same snapshots in 
time as in Figure [3] The actual equatorial internal energy cannot 
be precisely determined from the simulation because we only follow 
height-integrated quantities. 



significantly and compress the core. On the other hand, 
material outside of 1.3 M© spreads outwards consider- 
ably, so its internal energy drops via adiabatic expansion. 

While the main qualitative results of this simulation 
should be robust — the outwards transport of angular 
momentum, the conversion of rotational support into 
thermal pressure support, and the spreading of the outer 
regions — the quantitative results are only approximate. 
This is due to the height-integrated approximation, the 
confinement of mass to cylindrical shells, and the hard in- 
ner surface boundary condition. For example, as angular 
momentum is removed from the inner regions of the sim- 
ulation, some mass should flow towards the poles of the 
degenerate core at r cy \ = 0. This approach to spherical 



symmetry, and other multi-dimensional results, await a 
future simulation with the necessary physics (J. Schwab, 
et al., in preparation). 



4. THERMAL EVOLUTION: 10 3 - 10 4 



yr 



After the merger remnant has viscously evolved to 
a slowly rotating shear-free equilibrium, its subsequent 
evolution is determined by the heat generated previously 
during the dynamical and viscous phases of the merger. 
If the post-viscous temperature peak is sufficiently high 
(> 6 x 10 8 K at a density of 5 x 10 5 g cm -3 ), a con- 
vective C-burning shell will develop immediately follow- 
ing the viscous evolution on a timescale of ~ 10 6 s, 
prior to any significant thermal transport. While we 
are unable to extract precise temperatures from our one- 
dimensional viscous simulation, our preliminary multi- 
dimensional results suggest that this is indeed marginally 
the case for our fiducial 0.6 + 0.9 M© example. Such off- 
center convective burning converts the degenerate C/O 
core into a C-burning star from th e outside in, eventually 
yielding a degenerate O/Ne WD (fNomoto fc lb enl [19851 : 



iTimmes et all 119941: ISaio fc Nomotol 119981 ). During the 
~ 10 4 yr it takes for the burning layer to diffuse into 
the center, the shell-burning merger remnant radiates a 
small fraction of the nucl ear luminosity and ap pears as 
a near-Eddington source (jSaio fc Nomotolll998l ). 

For binary components of lower masses and thus lower 
gravitational potentials, the peak temperature following 
the viscous evolution will be lower, and convective C- 
burning will not immediately take place. In these cases, 
the further evolution will be determined by the transport 
of heat through the hot envelope. Since the envelope con- 
sists of material that was initially only halfway bound to 
the undisrupted WD, the conversion of its gravitational 
energy to heat and expansion work implies that the mate- 
rial is at least somewhat virialized, with a pressure scale 
height at the envelope's base that is a non-negligible frac- 
tion of the degenerate core's radius: h ~ i? CO ro- 

A simple order of magnitude estimate implies that the 
ratio of radiation to gas pressure in a virialized envelope 
is 



rad 



0.; 



Mr, 



Mr, 



M(7 

Men 



(G) 



where fi is the mean molecular weight. For double C/O 
WD mergers, this implies that radiation pressure is at 
least comparable to ideal gas pressure in the envelope^ 
Thus, the remnant radiates at a significant fraction of 
the Eddington luminosity while entropy is redistributed 
within the envelope. This occurs on a thermal timescale 
that is of order 



ithc 1' 



GA^qqj-q Ivl Gn v I -Rcorc 

4vrG(M corc + M cnv )c/ K 



10 4 yr 



0.2 cm 2 /g/ V R. 



10 cm 



(7) 



for core and envelope masses M COIC = 0.8 Af© and 
M C nv = 0.6 Mq, This timescale is equivalent to the 

7 However, this is not the case for lower mass cores and lower 
mean molecular weights (e.g., a 1 Mq star on the RGB or low mass 
double He WD merger remnants). 
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Figure 6. Snapshots of the T — p profile during the thermal evolu- 
tion of a 0.6 + 0.8 Mq merger remnant (thick solid lines). Profiles 
are equally spaced by 1000 yr and begin 500 yr after the end of 
the viscous evolution. Convective C-burning begins at 5790 yr at 
a mass coordinate of 0.836 Mq, which is represented by a square 
for each snapshot. The mass coordinate 0.735 Mq is represented 
by bullets. The dashed lines marks contours of constant specific 
entropy 10 80 and IO 8 4 erg g" 1 K" 1 as labeled. The thin solid 
lines denote P ra( j = Ptot/10 and Ptot/2. The dotted line marks the 
contour for which heating from C-burning equals neutrino cooling. 

Kelvin-Helmholtz contraction timescale for the material 
at the base of the envelope near R CO re, which transfers 
its entropy into regions further out in the envelope in the 
process of contracting onto the core. 

To explore this phase of thermal transport in more de- 
tail, we utiliz e the spherically sym metric stellar evolution 
code MESA (|Paxton et al. 20110. Because our height- 
integrated cylindrically symmetric viscous evolution (see 
Sj3]) cannot produce spherically symmetric results, we be- 
gin our thermal evolution with rough, but qualitatively 
correct, initial conditions: a degenerate 0.8 M© core at 
an initial temperature of 10 8 K surrounded by a high 
entropy 0.6 Mq envelope that does not exert significant 
pressure on the core. A primary of 0.8 Mq was chosen 
for its relatively lower gravitational potential, so that the 
maximum temperature in the post-viscous merger rem- 
nant would be lower than in the 0.9 + 0.6 Mq case. As 
previously mentioned, multi-dimensional viscous simu- 
lations are underway with ZEUS (J. Schwab, et al., in 
preparation); future work will use these results as ini- 
tial conditions for a more quantitatively correct thermal 
evolution. 

Given our poor understanding of mass loss in 
radiation-dominated H- and He-deficient envelopes, we 
turn off the options for mass loss in MESA. Addition- 
ally, the convective mixing length parameter is increased 
by a factor of 10 in the outer 0.1% of the remnant's mass 
in order to avoid severely pathological density profiles 
near the photosphere of the hot envelope, which oth- 
erwise force the code to take tiny timesteps. We will 
address these simplifications in more detail in 

Figure [6] shows T — p profiles (thick solid lines) at six 
snapshots in time equally spaced by 1000 yr, starting 

8 http://mesa.sourceforge.net/ (version 3635) 
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Figure 7. The total (solid lines) and radiative (dashed lines) 
luminosities vs. mass at the same snapshots in time as in Figure 
[6] The Eddington luminosity is 3.5 X 10 38 erg s -1 for a remnant 
mass of 1.4 Mq and k = 0.2 cm 2 g — 1 . 
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Figure 8. Temperature and radius vs. mass profiles at the same 
snapshots in time as in Figure [6] 

500 yr after the thermal evolution begins. The dotted 
line marks the contour for which heating from C-burning 
matches neutrino cooling. The thin lines show conditions 
for which the radiation pressure, P ra( j, equals one- half 
and one-tenth of the total pressure, Ptot- Bullets cor- 
respond to the trajectory of a fluid element at a mass 
coordinate of 0.735 M Q ; squares correspond to a mass co- 
ordinate of 0.836 Mq. The dashed lines demarcate con- 
tours of constant specific entropy at 10 80 and 10 8 4 erg 
g _1 K . As the squares show, fluid elements near the 
temperature peak lose entropy during the thermal evo- 
lution. A small amount of this entropy diffuses inwards 
into the degenerate core, reaching a mass coordinate of 
0.735 Mq after 6000 yr. However, since the diffusion 
timescale to the center of the degenerate core is ~ 10 6 
yr, the majority of the entropy from the inner regions is 
instead transferred outwards through the envelope. As 
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this material near the temperature peak loses entropy, it 
undergoes a Kelvin-Helmholtz contraction, falling from 
the hot extended envelope to the cold degenerate core 
and compressing to higher densities. This compression 
more than compensates for the loss of entropy and causes 
these fluid elements to increase in temperature. 

After 5790 yr have passed, the peak temperature 
reaches 6.3 x 10 8 K at a density of 4 x 10 5 g cm~ 3 , 
which is high enough for nuclear burning to overcome 
neutrino cooling and radiative diffusion, resulting in the 
onset of off-center convective C-burning at a mass coordi- 
nate of 0.836 Mq. Thus, even if convective shell-burning 
is avoided immediately following the viscous evolution 
and prior to significant thermal transport, it eventually 
sets in due to the cooling and compression of the base of 
the hot envelope. 

The key property of the post-viscous evolution that 
determines the time at which off-center convective burn- 
ing begins is the conditions at the temperature peak. 
We have chosen a somewhat arbitrary starting point for 
our thermal evolution calculation, with an initial tem- 
perature maximum of 6.0 x 10 8 K at 10 5 g cm -3 . If 
instead the post-viscous temperature peak were located 
at a higher density of 3 x 10 5 g cm -3 , only ~ 1000 yr 
would pass before convection begins. 

During the thermal evolution, the envelope approaches 
a constant luminosity profile at roughly one-third the Ed- 
dington limit, which is 3.5 x 10 38 erg s _1 for a remnant 
mass of 1.4 Mq and k — 0.2 cm 2 g _1 . In the interior of 
the envelope, most of this flux is carried by convection, so 
that the radiative luminosity profile is not constant (see 
the dashed lines in Fig. [7]). As entropy is redistributed 
throughout the envelope to yield a roughly constant total 
luminosity, shown as solid lines in Figure [3 the envelope 
spreads outwards significantly into a giant configuration 
until the majority of the envelope mass is > 10 11 cm. 
Mass-radius profiles are shown in the bottom panel of 
Figure [5] The top panel of Figure [5] shows that temper- 
atures throughout the extended regions of the envelope 
are < 10 6 K due to this expansion. 

Given the photospheric radius of 10 12 — 10 13 cm and 
the near-Eddington luminosity, the effective temperature 
during this cooling phase ranges from 5000 — 30000 K. 
However, we again emphasize that our artificial enhance- 
ment of the convective mixing length near the photo- 
sphere and the restriction on mass loss may change these 
results. While the merger remnant will radiate near the 
Eddington limit for ~ 10 4 yr, predictions of its tempera- 
ture and atmospheric properties will likely change with a 
more physical treatment of mass loss and convection. Fu- 
ture work on the thermal evolution of merger remnants 
will better quantify the appearance of these systems dur- 
ing the cooling and shell-burning phases. 

5. COMPARISON TO EARLIER WORK 

Two previous studies have examined the viscous and 
thermal evolution of the merger remnant using appropri- 
ate initial co nditions set by SP H simulations of double 
WD mergers. lYoon et al.l (l200l map the output of their 
SPH simulations to a one-dimensional stellar evolution 
code that includes the effects of rotation ()Yoon fc Langerl 
I2004D . Their initial remnant consists of the material pri- 
marily supported by pressure, which includes the mate- 
rial that had been shock-heated during the merger. The 



material primarily supported by rotation is assumed to 
reside in a Keplerian disk from which matter is accreted 
at a prescribed constant rate. Their viscous evolution, 
which occurs concurrently with the thermal evolution, 
includes the effects of hydrodynamic viscous transport 
processes but assumes that magnetic stresses are negligi- 
ble. Furthermore, the specific angular momentum of each 
shell is assumed to decrease on a parametrized timescale 
without liberating kinetic energy. 

Se veral impor t ant d ifferences distinguish our model 
from lYoon et al.l (|2007l ) 's. although some aspects are sim- 
ilar. The most significant physical difference is their re- 
striction to hydrodynamic angular momentum transport 
processes. If magnetic stresses are assumed to be neg- 
ligible, the timescale for angular momentum transport 
becomes comparable to the thermal timescale. As ma- 
terial from their rotating hot envelope cools and moves 
inwards, rotational support plays an increasingly impor- 
tant role and slows the contraction to the relatively slow 
angular momentum transport timescale, allowing more 
time for neutrino cooling and inwards thermal diffusion. 
Another important difference is their assumption that 
the material primarily supported by rotation following 
the merger resides in a Keplerian disk during the rem- 
nant evolution. We show that, instead, most of this ma- 
terial is converted into a high-entropy extended envelope 
during a phase of viscous evolution that occurs prior to 
the thermal evolution ( s ee Fig . [3]). 

Several of lYoon et all (|2007f )'s models assume accretion 
from the disk is negligible and that angular momentum 
loss from the total remnant occurs rapidly, so that an- 
gular momentum transport within the remnant is unim- 
portant (e.g., sequences Sa2, Aa2, and Ba2)0 Our model 
resembles these calculations if the total binary mass is as- 
sumed to be equal to their central remnant mass and if 
the heat generated b y viscous dissipati on is ignored. Due 
to these similarities. lYoon et al.l (|2007l ) also find that off- 
center C-ignition occurs in these models because of the 
relatively rapid Kelvin-Helmholtz contraction of the hot 
e nvelope surrounding t he de generate core. 

Ivan Kerkwiik et "all ((2010) consider the viscous evo- 
lution of mergers of equal-mass WDs in which 
both WDs are tida l ly di srupted, utilizing data from 
lLoren-Aguilar et al.l (|2009D 's SPH simulations. The re- 
sulting remnant has a temperature profile that peaks at 
the center, unlike remnants in which only one WD is dis- 
rupted, which have a temperature peak in material at 
the edge of the degenerate core. As in our work, their 
model assumes that magnetic stresses transport angular 
momentum in the outer materia l very rapidly. 

However, Ivan Kerkwiik et "aLI (2010) neglect the con- 
version of the rotational energy into heat, so that this 
material is rapidly transformed into degenerate material 
that exerts pressure on the core, yielding explosive C- 
burning at the center. While we agree that rapid vis- 
cous evolution occurs without significant heat transport 
between fluid elements, local conservation of energy de- 

9 Their test model Tal has no initial rotation nor accretion and 
would be very similar to our model, but the material at Tal's initial 
temperature peak is degenerate, with kT/Ep = 0.03, where Ep is 
the Fermi energy. Our model begins much more non-degenerately 
because of previous viscous heating, with kT/Ep = 1 at the tem- 
perature peak, and thus it experiences a much stronger compres- 
sional evolution. 
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mands that viscous dissipation converts the nearly Ke- 
plerian rotational energy into an equivalent amount of 
heat and expansion work. The resulting thermal pres- 
sure support prevents the outer material from compress- 
ing the central core until it undergoes Kelvin-Helmholtz 
contraction to rid itself of this excess heat. Furthermore, 
as mentioned in it is unclear that SPH simulations 
with accurate initial conditions allow for the tidal disrup- 
tion of both WDs, even if they have equal masses. Be- 
cause mass transfer in more accurate simulations begins 
at a larger orbital separation and persists for multiple or- 
bital periods prior to the merger, there is opportunity for 
symmetry breaking and the possibility that only one of 
the WDs is disrupted. Ongoing work will help to explore 
these outcomes in detail (M. Dan, et al., in preparation). 

6. CONCLUSIONS 

We have examined the post-merger evolution of un- 
equal mass double C/O WD systems taking into account 
the magnetic stresses that redistribute angular momen- 
tum rapidly compared to the thermal evolution of the 
remnant. This divides the merger evolution into a hier- 
archy of timescales: a rapid dynamical phase (10 2 — 10 3 
s) in which the less massive WD is tidally disrupted; 
an intermediate viscous evolution (10 — 10 8 s) during 
which the rotating envelope spreads and viscously heats 
without cooling; and, if convective C-burning does not 
set in first, a secular thermal phase (10 3 — 10 4 yr) in 
which the hot extended envelope radiates and contracts. 
This model is in contrast to earlier work that assumed 
the tidally disrupted WD is transformed into a Keplerian 
accretion disk fro m which the more massive WD accretes 
at a constant rate (iNomoto fe Ibenlll985l ; rSaio fe Nomotol 
[TflM [2?)0l IPiersanti et al.H2003allbn . " 

However, while the sequence of events is different, the 
outcome is qualitatively similar to previous work. In our 
models that do not immediately begin convective burning 
following the viscous evolution, thermal transport within 
and outwards from the envelope causes the base to be 
compressed, while in accretion models, accreting mate- 
rial is responsible for this compression. In both cases, 
however the rate of compression is set by the Eddington 
luminosity limit. Because of this relatively rapid com- 
pression and the prior dynamical and viscous heating, 
off-center convective C-burning is initiated and no SN la 
follows. Instead, the more likely outcome is the forma- 
tion of a high- mass O/Ne WD and possibly the collapse 
of the WD to a neutron star, as argued in previous calcu- 
lations. Although this is our preliminary conclusion, the 
physical model presented here differs significantly from 
that in previous work, so there are a number of unre- 
solved issues that need to be explored in more detail 
before the fate of WD mergers can be definitively estab- 
lished. 

In particular, the actual thermal evolution of the 
merger remnant may be strongly influenced by mass loss, 
which we have neglected in this study. Dust-driven, line- 
driven, and, given the importance of radiation pressure 
in the Eddington-limitcd envelope, continuum radiation- 
pressure-driven winds may play a significant role. In ad- 
dition to determining the thermal evolution of the inner 
layers, these winds may also influence the spectrum of a 
subsequent SN la. The escape velocity for material at 
the outer edge of the remnant (~ 10 13 cm) during the 



thermal evolution is ~ 60 km s~ 4 . If winds and/or the 
mass ejected during the coalescence and the viscous evo- 
lution propagate outwards for 10 4 yr, material will fill the 
surrounding volume out to ~ 2 x 10 18 cm at the time of 
an explosion. Such an outflow co uld possibly explain the 
variable fe.g.. iSimon et all 120091) or preferentially blue- 
shifted ([Sternberg et al.ll2011l ) Na I D absorption lines 
seen in some SNe la, although further detailed study of 
the mass loss mechanism is necessary to make quantita- 
tive predictions. 

Due to this neglect of mass loss, the eventual fate of 
the merger remnant remains uncertain. However, if the 
remnant collapses to a neutron star or yields a SN la, the 
resulting light curves will be affected by the presence of 
the extended envelope if it persists until the explosion. 
As ejecta from either of these events runs into ~ 0.1 Mq 
of material at radii of ~ 10 13 cm, the resulting shock ef- 
ficiently converts kinetic energy into thermal energy that 
will then be radiated on a photon diffusion timescale of 
several days. This makes a merger's collapse to a neu- 
tron star easier to observe; instead of reaching maximum 
light at ~ 1 d, such "enshrouded" explosions would peak 
at several days and decay from peak more slowly as well 
(jDarbha et alJl2010t IMetzger et al.ll2009bt ). However, in 
the case of SNe la, such a large amount of observable 
shock-deposited energy has been ruled out in at least 
one case by recent early-time measurem ents of SN 201 lfe 
(|Nugent et al.l[20lH iBloom et all 120121 ) . If SNe la arise 
from double WD mergers, their extended envelopes must 
not persist until the time of explosion. 

Future work is required to explore these outcomes, as 
well as to correctly calculate the envelope's structure 
near the photosphere. Our artificial enhancement of the 
convective mixing length near the surface allows us to 
overcome pathological density profiles, but it prevents 
us from reliably calculating the effective temperature of 
the remnant. Even so, we can make several robust pre- 
dictions from our model. Double C/O WD merger rem- 
nants will radiate near the Eddington limit for > 10 4 yr 
with bolometric luminosities ~ 10 38 erg s _1 . A defini- 
tive statement regarding their spectra and colors during 
their thermal evolution cannot be made without better 
determination of the mass loss rate. If outflowing dust 
does indeed play a significant role, they may be infrared- 
br ight and might b e identified in a sample such as that 
of iTisserandl (|201lD . Given the double WD merger rate 
of a f ew per 1000 yr (jNelemans et alj|2001t iRuiter et all 
2009), a Milky Way- type galaxy should harbor > a few 
dozen of these bright H- and He-deficient sources at any 
time. 

In this work, we have focused primarily on the merger 
of two C/O WDs. However, the same physical picture 
will also describe mergers for which one or both of the 
components is a He WD. Double He WD mergers will 
undergo a very similar evolution as we have described 
for C/O WD mergers, with appropriately adjusted en- 
ergy and timescales: a He-burning shell will ignite off- 
center and diffuse inwards, possibly with epis odic flashes 
due to the inefficienc y of neutrino cooling {iben 1990; 
iSaio fe Nomotol Il998t ). Once the burning shell reaches 
the center, the merger remnant appears as a core He- 
burning sdB or sdO star and eventually cools to become 
a C/O WD. The merger of a He and a C/O WD should 
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yield a He shell burning R CrB star ()Webbinkl 119841: 
iTbenl [19901 IClavtonl 1X9961) . although the final He shell 
flash of cooling post-AG B stars has also b een proposed 
as an R CrB progenitor (jlben et al.lll996[ ). Simulations 
of these mergers within the framework of our updated 
physical model will allow us to connect t heoretical pre- 
dictions with observational constraints ( Asplunc fet al.l 
[19971 IClavton et alJl2007L l20ll IJefferv et all 1201 lF and 
help to differentiate between these progenitor scenarios. 
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APPENDIX: VISCOUS SIMULATION DETAILS 

In this Appendix, we explain the computational details 
of our height-integrated viscous simulation whose results 
are shown in £j3j The approximations outlined there yield 
the Lagrangian time derivatives: 



~dt 



S d(r cy \v r ) 
dr cy \ 



fcyl 



dr, 



cyl 



dt 



(8) 
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where the cylindrically enclosed mass is m. 
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The height-integrated quantities £ and n are defined 
at zone centers, while the mass, radius, v r , and f2 are de- 
fined at zone boundaries. The radius at the inner bound- 
ary of the grid, where the enclosed mass is 0.9 M Q , is as- 
sumed to be fixed; i.e., the degenerate core is assumed to 
be incompressible. The angular velocity gradient across 
the inner and outer boundaries is set to zero. These 
boundary conditions, while clearly approximations, do 
conserve the total energy and angular momentum in the 
grid to better than 1% during the simulation. 

We approximate the radial gravitational acceleration 
a t (r cy i,z), where z is the height from the midplane, as 
the value at the midplane, f r (r cy \,z) ~ -Gm/r c 2 yl . This 
is equivalent to the assumption that pz 2 dz <C Sr^ yl . 

High-freque ncy spurious noise in v r is dampe d with 
the standard Ivon Neumann fc Richtmyerl (|1950l ) artifi- 
cial viscosity prescripti on modified for application t o 
curvilinear coordinates (|Tscharnuter fc Winkler! 11979ft . 
We define an artificial viscous pressure term as 



pl 2 V ■ v 




\5 i3 V ■ v 



for V • v < 



1,3) 



for V • v > 



and we add the terms —V • Q and — Q ■ (Vv) to the mo- 
mentum and entropy equations, respectively. The sym- 
metrized velocity gradient tensor is 



dxi 



(14) 



and Sij is the standard Kronecker delta,. The parameter 
I is a smoothing length, chosen to be the minimum zone 
spacing at each time step, and is constant throughou t 
the grid to preserve Tr[Q] = (jStone & Normanlll992T) . 



The condition on V • v in equation (|13l) ensures that the 
artificial viscosity only acts during compression. 
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